classdef Helper< handle
    properties (Constant)
        leaveOne = false %%要修改CalcAllErrors 
    end
    methods(Static = true)        
        %%X(204, 24), y(6,24)
        function showCorrelationEachBand(X, y, bandX)
            lenY = size(y, 1);
            
            figure
%             tiledlayout(2, 3);
            hold on
            
            cr = {'k', 'b', 'r'}
            for i = 1 : lenY
               yy = y(i, :);
               
               lenX = size(X, 1);
               corss = zeros(lenX, 1);
               
               for j = 1 : lenX
                  xx = X(j, :);
                  tmp = corrcoef(xx, yy);
                  corss(j) = tmp(1,2);
               end
%                nexttile
               [max(corss), min(corss)]
               plot(bandX, corss, cr{i});
            end            
            hold off
            legend('AN', 'AP', 'AK', 'location', 'east')
            xlim([bandX(1), bandX(end)]);
            xlabel('Wavelength(nm)')
            ylabel('Correlation coefficient')
        end
        
        %%画出每个样本的两个采样曲线
        function plotPair(XX)
            colors = [[1, 0, 0]; [0, 1, 0]; [0, 0, 1]; [0, 0, 0]; 
                      [1, 0, 1]; [0, 1, 1]; [0.5, 0.5, 0.5]; [0.7, 0.2, 0.2];
                      [0.3, 0.7, 0.1]; [0.3, 0.2, 0.7]; [0.7, 0.7, 0.1]; [0.7, 0.2, 0.7];
                      ];
                  
            lenX = size(XX, 2)/2;
            figure;
            hold on
            for i = 1 : lenX
                plot(XX(:, i), 'Color', colors(i, :));                
                plot(XX(:, i+lenX), 'Color', colors(i, :));
            end
            legend('1', '1', '2', '2', '3', '3', '4', '4', '5', '5', '6', '6')
            hold off
        end
        
        
        function plotTest()
            x=-3:.25:5;
            y=-0.3*x+3.5.*x.^2-x.^3+20*rand([1,length(x)]);
            [p,S]=polyfit(x,y,3);
            % 计算以p为系数的多项式在 x 中各点处的拟合值。将误差估计结构体指定为第三个输入，
            % 以便polyval 计算标准误差的估计值。标准误差估计值在 delta 中返回。
            
            [y_fit,delta]=polyval(p,x,S);
            % 绘制原始数据、线性拟合和 95% 预测区间 y±2Δ。
            uy=y_fit+2*delta;
            dy=y_fit-2*delta;
            % 绘制原始数据
            plot(x,y,'rx','LineWidth',1.2)
            hold on
            % 绘制拟合曲线
            plot(x,y_fit,'Color',[82,124,179]./255,'LineWidth',1.5)
            % 绘制置信区间
            plot([x',x'],[uy',dy'],'Color',[82,124,179]./255,'LineWidth',1.2,'LineStyle','--')
            fill([x,x(end:-1:1)],[uy,dy(end:-1:1)],[82,124,179]./255,'EdgeColor','none','FaceAlpha',.2)
            
            title('Linear Fit of Data with 95% Prediction Interval')
            legend('Data','Linear Fit','95% Prediction Interval')
            
        end
        function plotRegion(x, y, std)
            figure
            hold on
            plot(x, y, 'r');
            uy = y + std;
            dy = y - std;
            
            plot(x, uy, 'b--')
            plot(x, dy, 'b--')
            fill([x;x(end:-1:1)], [uy;dy(end:-1:1)],[82,124,179]./255,'EdgeColor','none','FaceAlpha',.2)
            
            hold off
            
        end
        
        %%% 画出12条光谱曲线
        function plot12Curve(xx, sTitle)
            figure
            hold on
            for i = 1 : length(xx)
                x = xx(i, :);
                plot(x)
            end
%             legend()
            title(sTitle)
            hold off
        end
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        function x_data =  extractMatrix(xx, comp)
            sz1 = size(xx{1},1);
            
            x_data = zeros(sz1, length(xx));
            for i = 1 : length(xx)
                x = xx{i};
                
                x_data(:, i) = x(:, comp);
            end
        end
        %%添加X轴标签
        function AddXTickLabels()
           xticklabels({'氮', '磷', '钾'})
        end
        
        function AddLegends()
            xticklabels({'氮', '磷', '钾'})
        end        
        
        function data = SaveStatsToFiles(fileSave, ret1, ret2)
            % % % 代码OK，但要增加MAPE，不需要R2total
            % % %             data = [arg.rmse; arg.R2; arg.rsd; arg.rpd];
            % % %             save(fileSave, '-ascii', '-tabs', 'data');
            % % %             R2total = arg.R2total;
            % % %             save(fileSave, '-ascii', '-tabs', '-append', 'R2total');
            
            data = [ret1.rmse; ret1.R2; ret1.rsd; ret1.rpd; ret1.mape];
            
            if nargin >2
                data2 = [ret2.rmse; ret2.R2; ret2.rsd; ret2.rpd; ret2.mape];
                data = [data; data2];
            end
            
            if isempty(fileSave) == false
                save(fileSave, '-ascii', '-tabs', 'data');
            end
        end
        
        function SaveValuesToFile(fileSave, ret1, ret2)
            data =  [ret1.y, ret1.yfit];
            save(fileSave, '-ascii', '-tabs', 'data');
            
            if nargin > 2
                data =  [ret2.y, ret2.yfit];
                fileID = fopen(fileSave, 'a');
                fprintf(fileID, '\r\n\r\n');
                fclose(fileID);                
                
                save(fileSave, '-ascii', '-tabs', 'data', '-append');
            end
        end
        
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        % % % % % % % % % % % % % % % % % % % % % % % % % % % %         %%
        function ret = CalcAllErrors(y, yfit, bShow, cc)
            if nargin < 4
                ret = Helper.CalcAllErrorsPatch(y, yfit, bShow);
                ret.y = y;
                ret.yfit = yfit;
                return
            end
            %%只针对 24折进行处理， 不是24折的，暂不处理
            if cc.NumTestSets ~= 24
               disp('error in 24-fold') 
            end
  
        end
        
        %% 想改但有没改成 2024-02-02
        function ret = CalcAllErrorsPatch(y, yfit, bShow)
             %% y =(24, 6)
             nn = size(y, 1);
             meanY = mean(y);
             rmse = rms(y - yfit);
             rmse_std = std(y - yfit);
             rsd = rmse ./ meanY;%%标准偏差系数 relative standard deviation
             %% rsd 应该给错了公式 2024-02-01， 应该是标准差/平均值
             
             
             %            rmse = rmse * sqrt(nn / (nn -1));
             
             %             R2total = norm(yfit - meanY)^2 / norm(y - meanY)^2;
             %             R2  = vecnorm(yfit - meanY).^2 ./ vecnorm(y - meanY).^2;
             R2total = 1 - norm(y - yfit)^2 / norm(y - meanY)^2;
             R2  = 1 - vecnorm(y - yfit).^2 ./ vecnorm(y - meanY).^2;
             
             ss = std(y); %% 分母为N-1，
             %             ss = std(y,1); %% 2024-01-30 改为分母N，
             rpd = ss ./ rmse; %%Relative Percentage Difference  residual prediction deviation
             
             
             err = abs(y - yfit);
             mape = mean(err ./ y);
             mape_std = std(err ./ y);
             
             if bShow
                 figure
                 hold on
                 xx = 1: size(y, 2);
                 plot(xx, rsd, '-rs')
                 plot(xx, R2, '-bo')
                 %                plot(rpd, '-kd')
                 hold off
                 %            legend('ph值', '有机质', '有效磷', '速效钾', '全氮', '砷', '镉', '铬', '铜', '镍', '铅')
                 xticks(xx)
                 Helper.AddXTickLabels()
                 title(R2total)
                 legend('相对标准偏差', 'R2')
             end
             ret.rmse = rmse;
             ret.R2 = R2;
             ret.R2total = R2total;
             ret.rsd = rsd;
             ret.rpd = rpd;
             ret.mape = mape;
             ret.mape_std = mape_std;
             ret.rmse_std = rmse_std;
        end
       
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        function PlotStemErrors(y, y_hat)
            residuals = y - y_hat;
            figure
            
            stem(residuals)
            xlabel('Observation');
            ylabel('Residual');
            legend
            
        end
        
        function PlotPLSRComps(X, y, nn)
            [XL,YL,XS,YS, beta, PCTVAR] = plsregress(X, y, nn);
            Helper.PlotPLSRvar(PCTVAR, true);
        end
        %%画出Y的累计方差
        function PlotPLSRvar(PCTVAR, bNewFig)
            nn = size(PCTVAR, 2);
            if bNewFig == true
                figure
            end
            hold on
            %             plot(1:nn,cumsum(100*PCTVAR(2,:)),'-bo');
            plot(1:nn,cumsum(100*PCTVAR(2,:)));
            xlabel('Number of PLS components');
            ylabel('Percent Variance Explained in y');
            
            %             plot(1:nn,cumsum(100*PCTVAR(1,:)),'-rs');            
            hold off
        end
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       
        function ret = isModelling(step)
            step4 = 8;
            step5 = 16;     step6 = 32;    step7 = 128;     step8 = 256;

           ret =  bitand(step, step4 | step5 | step6 | step7 | step8);
        end
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %         
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %         
% % % 2025-02-02 针对随机森林遍历全部预处理与模型% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
        function ret = DoAllMetalsFeaAndPre(arg)            
            x_proc = arg.x_proc;
            y_true = arg.y_true;
            y_train = arg.y_train;
            y_test = arg.y_test;
            trainIdx = arg.trainIdx;
            testIdx = arg.testIdx;
            seed = arg.seed;%%随机森林种子
            wave = arg.wave;%%特征波段数量
            nUseEachRemove = arg.nUseEachRemove;     
            bShow = arg.bShow;
            
            withNoPre = arg.withNoPre; %%0-不需要无预处理，1-需加上无预处理
            withNoFea = arg.withNoFea;
            pres = arg.pres + withNoPre; %%预处理个数
            feas = arg.feas + withNoFea; %%特征选择个数
            ret = cell(feas, pres);
            
            modelType = arg.modelType;
            %%预处理
            for preType = 1 : pres
                disp(['+++预处理步骤：', num2str(preType)]);
                rng(seed)%%设置随机种子
                %%先进行预处理
                if withNoPre == 0
                    x_data = PreHelper.DoOnePreProcessing(x_proc', preType, arg);
                else
                    x_data = PreHelper.DoOnePreProcessing(x_proc', preType - 1, arg);
                end
                
                %%预处理后的操作，重新得到训练集与测试集
                x_data = x_data';
                x_train = x_data(trainIdx, :);
                x_test = x_data(testIdx, :);
                
                if nUseEachRemove == 1                
                    y_test = y_true(testIdx, :);
                    [y_test, x_test] = WaveHelper.GetInsFromPos(y_test, x_test, arg);
                end
                
                %%特征选择
                if withNoFea == 1
                    disp('--------无特征选择步骤：');
                    rng(seed)%%设置随机种子
                    [ret1, ret2] = MyModel.DoOneModel(x_train, y_train, x_test, y_test, modelType, arg);
                    
                    
                    data = Helper.SaveStatsToFiles('', ret1, ret2);
                    ret{1, preType} = data;
                else
                    for feaType = 1 : feas
                        disp(['--------特征选择步骤：', num2str(feaType)]);

                        rng(seed)%%设置随机种子
                        kopt = FeaSelect.GetOptimalFeature(x_data, y_true, wave, feaType, 1:3, bShow);

                        [x_train2, x_test2] = FeaSelect.getXdataFromKopt(kopt, x_train, x_test, arg);

                        rng(seed)%%设置随机种子
                        [ret1, ret2] = MyModel.DoOneModel(x_train2, y_train, x_test2, y_test, modelType, arg);

                        data = Helper.SaveStatsToFiles('', ret1, ret2);   
                        ret{feaType, preType} = data;
                    end
                end
            end
        end
        %%把保存的mat文件保存在txt文件
        function SaveResultsFromPreFea(data, saveFile)
            fileID = fopen(saveFile, 'w+');           
            formatStr = repmat('%f\t%f\t%f\t\t', 1, size(data, 2));
            formatStr = [formatStr, '\n'];
            
            for pre = 1 : size(data, 1)
                outData = [];
                
                for fea = 1 : size(data, 2)
                    outData = [outData, data{pre, fea}];
                end
                
                fprintf(fileID, formatStr, outData');%要转置保存
                fprintf(fileID, '\n');
            end
            
            fclose(fileID);
        end
        
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
        function ret = DoAllMetalsModelFeaturePreProcess(x_data, y_data, arg)
               
            methods = arg.methods;
            ret = cell(size(methods, 1), 1);
            
            for metal = 1 : size(methods, 1)%%元素个数
                
                preType = methods(metal, 1);
                modelType = methods(metal, 2);
                
                if preType == 0 && modelType == 0
                    disp('preType == 0 & modelType == 0')
                    continue; 
                end
                ret{metal} = Helper.DoOneMetalWithFea(x_data, y_data, metal, preType, modelType, arg);              
            end
        end

        function ret = DoOneMetalWithFea(X, y, metal, preType, modelType, arg)
            %%X，24*204， y，24*6       
            bShow = arg.bShow;            
        
            x_data = PreHelper.DoOnePreProcessing(X, preType, false);          
            
            arg.metal = metal;
            arg.wave = 10;%%10个波段选择
            arg.bEach = false;            
            
            ret = cell(arg.feaTypes, 1);            
            
            % if modelType ~= 4 || metal ~= 3
            %     return;
            % end
            feaTypes = arg.feaTypes;

            for feaType = 1 : feaTypes %2024-07-27 %%专门针对 3-5的特征方法进行
                disp(['----fea：', num2str(feaType), '----model:', num2str(modelType), '----metal:', num2str(metal)]);
                arg.feaType = feaType;
                ret{feaType} = MyModel.DoOneModel(x_data', y', modelType, bShow, arg);
            end
        end
        %%2024-05-01 设置的每一个元素与指定settting，返回opt
        function ret = DoOneMetalWithOneSetting(X, y, metal, preType, feaType, arg)
            bShow = arg.bShow;
            
            x_data = PreHelper.DoOnePreProcessing(X, preType, false);
            arg.metal = metal;
            arg.wave = arg.wave;
            arg.bEach = false;
            arg.ncomp = 10;%%%%%%%%有问题要改
            
            disp(['----fea：', num2str(feaType), '----preType:', num2str(preType), '----metal:', num2str(metal)]);
            
            kopt = FeaSelect.GetOptimalFeature(x_data', y', arg.wave, feaType, metal, bShow);
            ret = sort(kopt);
% %             ret = kopt{metal};
% %             ret = sort(ret);
        end
        
        function data = SetResultForOneMetal(ret)
            len = length(ret);%%=6 金属个数
            data = cell(len, 1);
            for i = 1 : len
                metal = ret{i};
                if isempty(metal) == true
                    disp(['metal ', num2str(i), ' is empty'])
                    continue;
                end
                
                sz = length(metal);%%特征方法个数
                
                data{i}.rmse = zeros(1, sz);
                data{i}.R2 = zeros(1, sz);
                data{i}.rsd = zeros(1, sz);
                data{i}.rpd = zeros(1, sz);
                data{i}.mape = zeros(1, sz);
                data{i}.mape_std = zeros(1, sz);
                data{i}.rmse_std = zeros(1, sz);
                for j = 1 : sz
                    data{i}.rmse(j) = metal{j}.rmse(i);
                    data{i}.R2(j) = metal{j}.R2(i);
                    data{i}.rsd(j) = metal{j}.rsd(i);
                    data{i}.rpd(j) = metal{j}.rpd(i);
                    data{i}.mape(j) = metal{j}.mape(i);
                    data{i}.mape_std(j) = metal{j}.mape_std(i);
                    data{i}.rmse_std(j) = metal{j}.rmse_std(i);
                end
            end
        end
        % % % % % % % % % % % % % % % % % % % % % % % % % %
        function ret = DoOneModelWithAllSetting(x_proc, modelType, arg)             
            ret = cell(1, arg.preTypes);
            
            trainIdx = arg.trainIdx;
            testIdx = arg.testIdx;
            y_train = arg.y_train;
            y_test = arg.y_test;
%             x_train = arg.y_test;
            bWait = isfield(arg, 'wait');
            
            nUseEachRemove = arg.nUseEachRemove;     
            seed = arg.seed;%%随机森林种子            
            
            for preType = 1 : arg.preTypes
                rng(seed)%%设置随机种子
                 
                x_data = PreHelper.DoOnePreProcessing(x_proc', preType, arg);                
                x_data = x_data';
                x_train = x_data(trainIdx, :);
                x_test = x_data(testIdx, :);                
                
%                 FileSave = fullfile(arg.fileSavePath, 'bbb.txt');
%                 save(FileSave, '-ascii', '-tabs', 'x_train');
  
                disp(['----preType：', num2str(preType), '----modelType:', num2str(modelType)]);
                [ret1, ret2] = MyModel.DoOneModel(x_train, y_train, x_test, y_test, modelType, arg);   
                
                data = Helper.SaveStatsToFiles('', ret1, ret2);   
                ret{preType} = data;
%                 if bWait == true               
%                     pause(arg.wait) 
%                 end
            end            
        end
        
        
        
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        %% 在平均光谱线上 画出特征波段 2024-07-26
        function DrawSelectionOnMeanSpectral(x_data, fsss, offStart, offEnd)
            data = mean(x_data);

            lenX = length(data) + offStart + offEnd;
            for i = 1:length(fsss)
                xx = fsss{i}+offStart;
                yy = data(xx);
                
                figure
                plot(1:lenX, data);
                hold on
                plot(xx, yy, 'o');
                hold off
            end
            xlim([1, lenX]);
            
        end
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %         
        %% 用作合并两个输出文件 2024-07-29
        function ret = MergeRetFiles(file1, file2, fileSave, nCut)
            ret1 = load(file1, 'ret');
            ret2 = load(file2, 'ret');
            ret1 = ret1.ret;
            ret2 = ret2.ret;
            
            if length(ret1) ~= length(ret2)
                disp('error in length')
                return
            end            
            ret = cell(length(ret1), 1);
            for i = 1 : length(ret1)
                data1 = ret1{i};
                data2 = ret2{i};
                ret{i} = [data1(1:nCut-1); data2(nCut : end)];
            end
            
            save(fileSave, 'ret');
        end
        
% % % % % % % % % % % % % % % %        
%%%%%%%% 绘制残差图与45度线
        function DrawResidua(y_true, y_fit, bShowPoint)
            if nargin <=2
                bShowPoint = false;
            end
                
            sizeY = size(y_true, 2);
            
            sTitle = {'N', 'P', 'K'};
            if bShowPoint == true
                len = size(y_fit, 1);
                txt = cell(len, 1);
                
                for kk = 1 : len
                    txt{kk} = num2str(kk);
                end
            end
            
            figure;
            tt= tiledlayout(1, 3);
            tt.TileSpacing = 'compact';
            tt.Padding = 'compact';
            for i = 1 : sizeY
% %                 figure,
                nexttile
                plot(y_true(:, i), y_fit(:, i), '.', 'markersize', 17);
                title(sTitle{i})
                a = max(max(y_fit(:, i)), max(y_true(:, i))) * 1.05;
                xlim([0, a])
                ylim([0, a])
                hold on
                plot([0, a], [0, a], '-')
                
                if bShowPoint == true
% %                     text(y_fit(:, i), y_true(:, i), txt);
                    text(y_true(:, i), y_fit(:, i), txt);
                end
                
                axis square
                ylabel('预测值');
                xlabel('真实值');
                hold off
            end
        end
        
        %%画出含量分布直框图
        function DrawHistgram(y_true)
            figure
            tiledlayout(1, 3);
            nexttile
            histogram(y_true(:, 1));
            title('N');
            
            nexttile
            histogram(y_true(:, 2));
            title('P');
            
            nexttile
            histogram(y_true(:, 3));
            title('K');
        end
        
        %%%合并样本针对 flag=2进行合并
        function [yy2, xx2] = ProcessOriginalData(y_data, x_data, flag, seq, index)
            ins = flag == 2;          
            
            y_data = y_data(ins, :);
            x_data = x_data(ins, :);
            flag = flag(ins, :);
            seq = seq(ins, :);
            index = index(ins, :);
            
            ind = split(index, '-');
            
            first = str2double(ind(:, 1));
            second = str2double(ind(:, 2));
            
            keys = unique(first);
            xx2 = [];
            yy2 = [];
            for i = 1 : length(keys)
                key = keys(i);
                
                inn = first == key;
                
                mmx = mean(x_data(inn, :));
                mmy = mean(y_data(inn, :));
                xx2 = [xx2; mmx];
                yy2 = [yy2; mmy];
            end
            
        end
        
        %%合并样本,正负样本合并成1个，
        function [yy2, xx2, flag, seq, index] = ProcessOriginalData2(y_data, x_data, flag, seq, index, merge)
            ins = flag == 1;          
            
            y_data = y_data(ins, :);
            x_data = x_data(ins, :);
            flag = flag(ins, :);
            seq = seq(ins, :);
            index = index(ins, :);
            if merge == false
                yy2 = y_data;
                xx2 = x_data;
                return
            end            
            
            ind = split(index, '-');
            
            first = str2double(ind(:, 1));
            second = str2double(ind(:, 2));
            
            keys = unique(first);
            xx2 = [];
            yy2 = [];
            for i = 1 : length(keys)                
                key = keys(i);
                
                inn = first == key;
                
                xx_data = x_data(inn, :);
                yy_data = y_data(inn, :);
                second2 = second(inn, :);
                
                key2 = unique(second2);
                for j = 1 : length(key2)
                    kkk = key2(j);
                    in2 = second2 == kkk;
                    if sum(in2)>1
                        mmx = mean(xx_data(in2, :));
                        mmy = mean(yy_data(in2, :));
                        xx2 = [xx2; mmx];
                        yy2 = [yy2; mmy];
                    end
                end
            end            
        end
        
        %% 保存误差，并按最值大小排序。
        function SaveSortDiffValue(fileSave, ret1, ret2)
            
            diff1 = abs(ret1.y - ret1.yfit);
            diff2 = abs(ret2.y - ret2.yfit);
            
            [vv1, ii1] = sort(diff1, 'descend');            
            [vv2, ii2] = sort(diff2, 'descend');
            
            save(fileSave, '-ascii', '-tabs', 'vv1');
            
            fileID = fopen(fileSave, 'a');
            fprintf(fileID, '\r\n\r\n');
            fclose(fileID);                

            save(fileSave, '-ascii', '-tabs', 'vv2', '-append');            
            
        end
        
        %%2024-11-02 找到最大的残差点
        function [savePos, saveM] = FindMaxDiff(y_true, y_fit, args)                
            klen = args.klen;        
            llen = size(y_true, 2);
            
            savePos = zeros(klen, llen);
            saveM = zeros(klen, llen);
            RR2 = zeros(klen, llen);
            for k = 1 : llen
                y = y_true(:, k);
                ypred = y_fit(:, k);                
                
               
                pos = true(size(y));
                ind = 1 : length(y);
                
                for i = 1 : klen
                    y = y(pos);  
                    ind = ind(pos);

                    ymean = mean(y);
%                     ypred = ypred(pos, :);
                    ypred = ypred(pos);

                    R2 = Helper.calcR2(y, ypred);
                    RR2(i, k) = R2;
                    
                    
% % % % %      OK：方法3,增加动态调整的均值              
                    nnn = length(y);
                    SSR = sum((y - ymean).^2) * ones(size(y));
                    SSE = sum((ypred - y).^2) * ones(size(y));
                    
                    ss = (sum(y) - y) / (nnn - 1);
                    dd = y.^2 - nnn * ymean.^2 + (nnn - 1) * ss.^2;
                    eachSSR = SSR - dd;
                    eachSSE = SSE - (ypred - y).^2;
                    ratio = eachSSE ./ eachSSR;
                    [M, in] = min(ratio);                    
                    
% % % % %      OK：方法2，增加删减后的 R2变化  
% % % %                     SSR = sum((y - ymean).^2) * ones(size(y));
% % % %                     SSE = sum((ypred - y).^2) * ones(size(y));                    
% % % %                     eachSSR = SSR - (y - ymean).^2;
% % % %                     eachSSE = SSE - (ypred - y).^2;
% % % %                     ratio = eachSSE ./ eachSSR;
% % % %                     [M, in] = min(ratio);
                    
% % % %                    OK,最早的方法1：
% % % %                     eachMean = (y- ymean).^2 + 0.000001;
% % % %                     eachSSE =  (ypred - y) .^2;
% % % % 
% % % %                     ratio = eachSSE ./ eachMean;
% % % %                     [M, in] = max(ratio);

                    pos = true(size(y));
                    pos(in) = false;          

                    savePos(i, k) = ind(in);
                    saveM(i, k) = M;
                end
            end         
            Helper.DrawDiffRR2(RR2)
            aa = [savePos, saveM];
            save(args.fileSaveMaxPos, '-ascii', '-tabs', 'aa');
            
            diff = RR2(2:end, :) - RR2(1:end-1, :);
            [B, I] = maxk(diff, 15);
            save(args.fileSaveDiffPos, '-ascii', '-tabs', 'I');
            
        end
        %%% 画出显著差异
        function DrawDiffRR2(RR2)        
            figure
            tiledlayout(1, 3);
            llen = size(RR2, 2);
            for k = 1 : llen
                 aa = RR2(2:end, k) - RR2(1:end-1, k);
                 
                 nexttile
                 plot(aa);
                 [B,I] = maxk(aa, 10)
            end
            title('DrawDiffRR2');
            
            figure
            tiledlayout(1, 3);
            llen = size(RR2, 2);
            for k = 1 : llen
                 nexttile
                 plot(RR2(:, k));                 
            end
            title('RR2');
        end
        
        %%计算R2
        function R2 = calcR2(y, ypred)
            ymean = mean(y);
            SSR = sum((y- ymean).^2);
            
            err = ypred - y;
            SSE = sum(err.^2);
            
            R2 = 1 - SSE / SSR;
        end
        %% 去除指定删除点
        function varargout = removePoints(savePos, varargin)
            szData = length(varargin);
            
            for k = 1 : szData
                y = varargin{k};
            
                zeroPos = true(size(y,1), 1);
                zeroPos(savePos) = false;

                varargout{k} = y(zeroPos, :);
            end  
        end
        
        
        
        %% 去除指定删除点，任意入参
        function varargout = removeYandX(savePos, varargin)
            szData = length(varargin);
            
            for k = 1 : szData
                y = varargin{k};
            
                zeroPos = true(size(y,1), 1);
                zeroPos(savePos) = false;

                varargout{k} = y(zeroPos, :);
            end  
        end

         function varargout = removeSamePoints(args, varargin)
            savePos = load(args.fileSaveMaxPos, '-ascii');
            
            nRemove = args.nRemove;
            
            pos = savePos(1:nRemove, 1:3);
            pos = unique(pos);
            
            szData = length(varargin);
            
            for k = 1 : szData
                y = varargin{k};
                
                zeroPos = true(size(y,1), 1);
                zeroPos(pos) = false;
                
                varargout{k} = y(zeroPos, :);
            end
            varargout{szData+1} = pos;
            
         end

         function ret = StatConcentration(y)           
             ret = zeros(5,3);
             ret(1, :) = min(y);
             ret(2, :) = max(y);
             ret(3, :) = mean(y);
             ret(4, :) = std(y);
             ret(5, :) = std(y) ./ mean(y);
             ret = ret';
             save('StatAvailables.txt', '-ascii', '-tabs', 'ret');
         end
         
         function ret = loadExcelAllData(dataSrcFile)
            data = readtable(dataSrcFile, 'sheet', '总表');

            ret = table2array(data(2:374, 4:227));
            
            p = randperm(size(ret, 1));
            a = p(1:223);
            
            ret = ret(a, :);
         end

             
    end
end